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Abstract 

We study by kinetic Monte Carlo simulations the dynamic behavior of a Ziff-Gulari-Barshad 
model with CO desorption for the reaction CO + O — > CO2 on a catalytic surface. Finite-size 
scaling analysis of the fluctuations and the fourth-order order-parameter cumulant show that below 
a critical CO desorption rate, the model exhibits a nonequilibrium first-order phase transition 
between low and high CO coverage phases. We calculate several points on the coexistence curve. We 
also measure the metastable lifetimes associated with the transition from the low CO coverage phase 
to the high CO coverage phase, and vice versa. Our results indicate that the transition process 
follows a mechanism very similar to the decay of metastable phases associated with equilibrium 
first-order phase transitions and can be described by the classic Kolmogorov- Johnson-Mehl-Avrami 
theory of phase transformation by nucleation and growth. In the present case, the desorption 
parameter plays the role of temperature, and the distance to the coexistence curve plays the role 
of an external field or supersaturation. We identify two distinct regimes, depending on whether 
the system is far from or close to the coexistence curve, in which the statistical properties and the 
system-size dependence of the lifetimes are different, corresponding to multidroplet or single-droplet 
decay, respectively. The crossover between the two regimes approaches the coexistence curve 
logarithmically with system size, analogous to the behavior of the crossover between multidroplet 
and single-droplet metastable decay near an equilibrium first-order phase transition. 

PACS numbers: 82.65.+r, 64.60.Ht, 82.20.Wt, 05.40.-a 
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I. INTRODUCTION 



The study of phase transitions and critical phenomena in nonequilibrium systems is a 
subject of great interest. In particular, the study of surface reaction models has attracted 
considerable attention 1]. These models not only exhibit rich and complex behavior, but 
they can also explain a wide range of experimental results associated with catalysis and 
could be very useful for designing more efficient processes. The potential applications of 
improved catalytic reactions are a powerful reason to pursue this line of research |2j . Unlike 
the decay of metastable phases near first-order equilibrium phase transitions, the decay 
of metastable phases near nonequilibrium phase transitions still lacks a well-established 
theoretical framework. 

The Ziff, Gulari, and Barshad (ZGB) model is a lattice-gas adsorption-reaction model 
that describes some kinetic aspects of the catalytic oxidation of carbon monoxide on a crystal 
surface The ZGB model assumes that the reaction between CO and O2 on a surface 
proceeds according to the Langmuir-Hinshelwood process: 

CO(g) + S -> CO (a) 
2 + 2S -> 20(a) 
CO(a) + 0(a) -> C0 2 (g) + 2S, 

where S is an empty site on the surface, and (g) and (a) refer to the gas and adsorbed 
phase, respectively. The process is controlled by a single parameter y, which represents the 
probability that the next molecule arriving at the surface is CO, i.e., it is proportional to 
the partial pressure of CO. The model exhibits two kinetic phase transitions, a continuous 
one at y = y 1 and a discontinuous one at y = y 2 , where yi < y 2 . When y < y x , all the 
sites become occupied by oxygen, the so-called oxygen-poisoned state. If y > y 2 , all the 
sites become occupied by CO molecules, the so-called CO-poisoned state. Real systems 
do not possess an oxygen-poisoned state because oxygen does not impede the adsorption 
of CO. However, transitions between states of low and high CO coverage 8qo (where 9^ 
is the fraction of surface sites occupied by CO) have been observed experimentally jj 
At low temperatures, as y increases, there is a discontinuous increase in #co; accompanied 
by a discontinuous drop in the C0 2 production rate. Above a critical temperature the 
discontinuities disappear, and the C0 2 production decreases continuously. This type of 
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behavior can be reproduced by modifying the ZGB model to include a CO desorption rate, 
, a model we for brevity will call the ZGB-k model [l(|. For this model, there 
is a distinction between high and low CO-coverage phases only for k below a critical value 
k c , while above k c the CO coverage varies smoothly with y. Thus, the transition value 1/2 
becomes a function of fc, corresponding to a coexistence curve yi {k) that terminates at the 
critical point y2{k c ) 111 1X21] . The ZGB-k model does not have a totally poisoned CO state, 



totally 

HQ. 



and hysteresis is observed in #00 as y is varied close to 2/2 (&) [id, This hysteresis is 
associated with well-defined metastable phases of the model. 

The main aim of this paper is to understand the dynamical response of the model near the 
discontinuous transition. We present an extensive finite-size scaling analysis of results from 
kinetic Monte Carlo simulations that indicates, more conclusively than previous studies, 
that the system undergoes a first-order nonequilibrium phase transition along a coexistence 
curve that terminates at a critical point. We calculate several points on the coexistence 
curve. Next we measure the lifetimes of the metastable phases associated with the decay 
from the low (high) CO coverage phase to the high (low) coverage phase. We find that the 
statistics of the metastable lifetimes are well described by the Kolmogorov-Johnson-Mehl- 



Avrami (KJMA) IjJ, [14], [l5( theory of phase transformation by nucleation and growth near 



a first-order equilibrium phase transition. 

The outline of the paper is as follows. In Sec. [H] we define the model and describe the 
Monte Carlo simulation techniques used. In Sec. IIHI we present and discuss the numerical 
results obtained: in Sec. lIII Al we show how we calculate the coexistence curve and present a 
finite-size scaling analysis of the fluctuations and of the fourth-order cumulant of the order 
parameter; in Sec. lIIIBl we present the measurements of the lifetimes of the metastable states 
associated with the transition and show how their behavior is described by the KJMA theory. 
Our conclusions are summarized in Sec. IIVI 



II. MODEL AND SIMULATION 



The ZGB model with desorption is simulated on a square lattice of linear size L that 
represents the catalytic surface. A Monte Carlo simulation generates a sequence of trials: 
CO or O2 adsorption with probability 1 — k and CO desorption with probability k. In the 
case of adsorption a CO or O2 molecule is selected with probability y and 1 — y respectively 
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FIG. 1: Schematic representation of the algorithm. See discussion in the text. 



|, These probabilities are the relative impingement rates of both molecules and are 
proportional to their partial pressures. The algorithm works in the following way. A site i is 
selected at random. In the case of desorption, if i is occupied by CO the site is vacated and 
the trial ends, if not the trial also ends. In the case of adsorption, if a CO molecule is selected 
it can be adsorbed at the empty site i if none of its nearest neighbors are occupied by an O 
atom. Otherwise, one of the occupied O neighbors is selected at random and removed from 
the surface, leaving i and the selected neighbor vacant. This move simulates the CO + O 
— > CO2 surface reaction following the adsorption of CO. O2 molecules can be adsorbed only 
if a pair of nearest-neighbor sites are vacant. If the adsorbed molecule is selected to be O2 
a nearest neighbor of i, j, is selected at random, and if it is occupied the trial ends. If both 
i and j are empty, the trial proceeds, and the O2 molecule is adsorbed and dissociates into 
two O atoms. If none of the remaining neighbors of i is occupied by a CO molecule, one O 
atom is located at i, and if none of the neighbors of j is occupied by a CO molecule, then 
the other O is located at j. If any neighbors of i are occupied by a CO, then one is selected 
at random to react with the O at i such that both sites are vacated. The same reaction 
happens at site j if any of its neighbors are filled with a CO molecule. This process mimics 
the CO + O — > CO2 surface reaction following O2 adsorption. A schematic representation 
of this algorithm is shown in Fig. ^ We emphasize that the ZGB model, both with and 
without CO desorption, is an intrinsically nonequilibrium model that is fully defined by 
these dynamic rules. In contrast to systems considered in equilibrium thermodynamics, its 
properties are not derived from a Hamiltonian. We shall return to this point in Sec. IIIIBI 

For our simulations we assume periodic boundary conditions. The time unit is one 
Monte Carlo Step per Site (MCSS), in which each site, on average, is visited once. For 
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FIG. 2: Contour plot of the projection of 10 MCSS of simulation onto the ternary phase diagram 
for k = 0.02, y = 0.5332, and L = 100. 

measurements of stationary quantities, the system was allowed to reach stationarity before 
data were recorded for analysis. Averages were taken over 10 3 independent simulation runs. 

III. RESULTS 

We use a standard ternary phase diagram to plot the fraction of sites occupied by CO 
molecules: the CO coverage, 6^0; the O coverage, 6*0, and the fraction of empty sites, #e- 
In Fig. |21 we present a contour plot of a histogram based on the projection of 10 6 MCSS 
onto the plane of the phase diagram. For the chosen parameters and observation times the 
system undergoes several transitions between the low and high CO coverage phases. From 
the fact that the set of phase points is nearly parallel to the line #o/#e = 1/2, it is evident 
that the CO coverage gives more information about the kinetic phase transitions than 60 or 

0E- 
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FIG. 3: Order-parameter probability distribution, P(8co), shown vs y for k = 0.03 and L = 100. 
The distribution for the value of y closest to the coexistence value is shown with a bold line. 

A. Determination of the Coexistence Curve 

We estimate P(9qo), the probability distribution for Oco, by recording the number of 
times Ni that the coverage fell in the intervals [0, e), [e, 2e), [1 — e, 1] (e = 0.01), such that 

Ni = N is the total number of MCSS. Then, the probability that 8qo has a value in the 
interval [(i — l)e, ie) is Pj = Ni/Ne, such that 



/ P{d co )dd co = Y j P l e = l. 
Jo 



(1) 



In Fig.Elwe show P($co) versus y. In the regions (below and above 1/2) where the histograms 
are unimodal, the system consists of one single phase. For a very narrow range of y, the 
histograms are bimodal, indicatin g tw o distinct phases. At the coexistence point y%{h), the 
areas under both peaks are equal 16. "v\. 

We define a measure of the fluctuations in #co in a L x L system in the standard way as 

Xl = L 2 ((61 o ) l -(0 co )1), (2) 

where 



? co)i 



9 n co P{e C o)d6 C o ■ (3) 
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FIG. 4: The order-parameter fluctuation measure Xl, shown vs y for A; = 0.02 and four system 
sizes. The dotted lines are guides to the eye. The values of Xl have an error of approximately 5%. 

We measure Xl as a function of y for a fixed value of k and several values of L. At a 
first-order equilibrium phase transition, the order-parameter fluctuations increase with the 
system size, such that the maximum value of Xl ~ L d , where d is the spatial dimension of the 



system 17LI1S. 
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201 ]. We will take the same scaling behavior to indicate a nonequilibrium 



first-order transition. 

In Fig. |3] we show Xl vs y for four system sizes at k — 0.02. For the four values of L 
used, Xl displays a clear peak, which moves and increases in height with increasing L. 
In Fig. E2a) we plot ln(X™ ax ) versus ln(L) for several values of k. A linear fit indicates a 



power- law divergence with L, such that the maximum value scales as X\ 



L d ' with d! 



1.96 ±0.02, 1.91 ±0.05, and 1.58 ±0.04 for k = 0.02,0.03, and 0.04, respectively. A different 
method to extract the power-law exponent, which has some advantage in eliminating the 
effects of a nonsingular background term (as in Xl = f + gL d ' with / and g constants), is 
to consider 



In 



ymax 
Ji£"max 



/In 6 = d' + 0(l/lnb) (4) 

with L fixed at a relatively small value (here, L = 20), and b > 1. For large L and b, 
the correction term is proportional to //(pin 6), so that the exponent can be estimated 
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FIG. 5: (a) Plot of ln^j? 1 **) vs ln(L) and (b) plot of \n(X™*/Xf^)/ ln(b) with L = 20 vs 
l/ln(6), both for several values of and including all four system sizes. X™ ax is the maximum 
value of Xl, taken from figures similar to Fig. The straight lines are the best linear fits to 
the data and give X™ ax ~ L d ' with (a) d' = 1.96 ± 0.02, 1.91 ± 0.05, and 1.58 ± 0.04; and (b) 
d! = 2.008 ± 0.005, 2.006 ± 0.008, and 1.52 ± 0.04; for k = 0.02, 0.03, and 0.04, respectively. 

by plotting the left-hand-side of Eq. (J3J) vs l/ln& and extrapolating to 1/ In b = 0, as in 
Fig.EKb). The resulting estimates are d' = 2.008 ± 0.005, 2.006 ± 0.008, and 1.52 ±0.04 for 
k = 0.02, 0.03, and 0.04, respectively. 

These results indicate that the system undergoes a first-order phase transition, d' ~ 
d = 2, at y = V2{k), generating a coexistence curve that terminates at a critical point 
0.03 < k c < 0.04. A previous estimate based on a study of the fractal dimension of the 
interface between the two phases gives 0.039 < k c < 0.04 9j. Another study, which estimates 
k c (L) as the value of k where the double-peaked histograms become single-peaked, gives 
k c (L — > oo) 



0.04060 



11 1 . However, from the known two-peaked shape of the order- 
parameter distribution at the equilibrium Ising critical point , we believe that this method 
should yield a slight overestimate of k c . Reference llj also reports preliminary results on 
the fourth-order cumulant of the CO coverage that are consistent with an Ising-like critical 
point at k c = 0.0406. However, the cumulants were calculated only for very small lattice 
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sizes (L = 10, 20, and 40) and do not constitute a definite proof of the location of the critical 
point, as the authors duly point out nj. Below we present a similar cumulant study, but 
with larger lattice sizes. 

We also calculated the relation between the value of the probability distribution P(6qo) 
at either of the peaks of the bimodal distribution, P max , and its value at the local minimum 
between the peaks, P mm . For a first-order equilibrium phase transition these quantities 
satisfy the relation 

-Pmax 

— — oc exp(cL) , (5) 

i min 

where c is proportional to the equilibrium interface tension between the two phases. If the 
relation is applied to the present system, we would expect c to be positive and decrease with 
increasing k. Our results indicate that c(k) is positive only for k < 0.03, suggesting that 
0.02 < k c < 0.03. These results corroborate again that the system has a first-order phase 
transition for small k. However, they suggest a much lower value for k c than indicated by 
our other techniques and the previous results by others 0, [ll|. We find this result quite 
interesting and believe it may be due to several reasons. Most obvious is the significant 
difficulty in locating and measuring the peaks, which are extremely narrow in y. Perhaps 
more significant is the fact that this non-Hamiltonian nonequilibrium system does not possess 
a well-defined surface tension that could be associated with the parameter c in Eq. (|SJ). This 
point will be further discussed in Sec. IIIIB1 where we also show that the cluster interfaces 
in this system are much rougher than in conventional Hamiltonian systems. As a result, we 
do not consider the method for determining k c , based on Eq. (j^J), very accurate. The fact 
that different techniques give different results are an indication of the difficulties associated 
with locating k c in this model, even for relatively large systems. 

A useful tool for detecting phase transitions in simulations of finite equilibrium systems 



is the fourth-order reduced cumulant of the order parameter 
the form 
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For 9 co it takes 



u L = l- —2, (6) 



3u 



/in = <(0CO - (#CO» n > = / (#CO " ^CO>) n P(#Co)^CO, (7) 



where, 



is the nth central moment of the CO coverage. The equal-area bimodal distribution corre- 
sponding to coexistence yields a positive maximum for the cumulant vs y, flanked on either 
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FIG. 6: Dependence of the fourth-order reduced cumulant ul on y, for (a) k = 0.02 and (b) 
k = 0.04. The P(9qo) histograms indicate that the minima of ul correspond to the transitions 
from unimodal to bimodal distribution. The maximum between them gives the coexistence point 
y 2 (k,L). The horizontal lines in the insets correspond to ul = 2/3 (dashed) and ^j sing ~ 0.61 
(dotted). 

side by negative minima and approaching zero far away from the transition. Since the cumu- 
lant essentially is a tool to determine the shape of the order-parameter distribution, it can 
also be used for nonequilibrium phase transitions, such as the one studied here. The maxima 
of ul define the L-dependent coexistence line, y 2 (k, L). In Fig.Elwe show the dependence of 
ul on y for several values of L and for k = 0.02 and k = 0.04, respectively. When k = 0.02, 
the maximum value of ul is very close to 2/3, however when k = 0.04 the maximum is very 
close to 0.61, consistent with the proximity of an Ising-like critical point. 

The finite-size scaling theory of equilibrium first-order phase transitions implies that the 
shift in the position of the transition in a finite system of linear size L with periodic boundary 
conditions is inversely proportional to the system volume, L d jisl. I20I I22 1 (here the dimension 
d — 2). Although there is no analogous scaling theory for the present nonequilibrium system, 
we here attempt to use the same scaling relation, 

y 2 (k,L)-y 2 (k)cxL- 2 , (8) 

where y 2 {k) is the transition value of the CO adsorption rate y in the infmite-L limit. In 
Fig. [3 we plot y 2 (L) vs 1/L 2 for L = 20,40,60,80,100,200, and 300 for several values 
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FIG. 7: Critical CO adsorption rate y2(k,L), shown vs 1/L 2 for several values of the desorption 
rate k. The straight lines are weighted least-squares fits, yielding the estimated y 2 (k) for L = oo. 



of k. The error bars are calculated as the half width of the peaks of the ul vs y curves 



119, 



201 ]. As seen from the figure, the points fall very close to the straight line representing 
a weighted least-squares fit, yielding a good estimate (within 1CT 2 % ) of y 2 {k) for each k 
when L — > oo. Our result y 2 = 0.5257(3) for k = is consistent with previous studies that 
give y 2 = 0.52560(1) [TJ and y 2 = 0.52583(9) [23] for L -> oo. 

In Fig. |H1 we show several points on the coexistence curve between the low and high CO 
coverage phases. The coexistence curve is almost linear with only a slight curvature. By 
extrapolating a quadratic fit of the data to k = 0.04060 we obtain y 2 = 0.542(3) in agreement 
with the previous result y 2 = 0.54212(10) for k c jll|. 



B. Metastability 

In this section we calculate the time associated with the decay from the low (high) CO 
coverage phase to the high (low) CO coverage phase, r p (r d ), and determine its dependence 
on the CO pressure and the system size. In order to do this, we prepare the system with 
initial pressure y = y w such that y\ < y w < y 2 {k) (y w > y 2 {k)) and then suddenly change 
y, such that y > y 2 {k) (y < y 2 {k)). Then the initial low (high)- coverage phase becomes 
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FIG. 8: Some points of the coexistence curve, analogous to the pressure vs temperature phase 
diagram for a fluid in equilibrium. The continuous line represents a quadratic fit to the data. The 
vertical dashed line indicates the estimate k c ~ 0.039 from Ref. fj]. 

metastable and eventually decays to the high (low)-coverage phase. The system is considered 
to have left the metastable phase once its coverage reaches a certain cutoff value 8q Q . To 
avoid that recrossing events back to the metastable phase are mistaken for decay events, 
the cut-off is selected such that it is not too close to the metastable coverage value. The 
statistics of the decay times are analyzed for n = 500 independent runs. 

Since the value of y w that determines how far the initial system is from the transition 
point 2/2 (&) is somewhat arbitrary, it is necessary to evaluate how the decay times depend on 
it. Figure |Hf a) indicates that, in the region of interest, the average decay time from the low 
to the high CO coverage phase, (r p ), while being dependent on the final pressure y, is fairly 
independent of the pressure at which the system is prepared, y w . In the following we then 
take y w = 0.45 to calculate r p . It is also necessary to evaluate the dependence of (r p ) on 
the selected cut-off value 0q O . In Fig. EJb) we plot (r p ) vs 8q for different values of k and 
y. Clearly, (r p ) increases with 0q O , however there is a region where it is relatively weakly 
dependent on the cutoff. We choose to perform our measurements of (r p ) at 0q O = 0.65, 
well inside this region. 

Figure HH indicates that the average decay times associated with the decontamination 
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FIG. 9: (r p ) as a function of (a) y w with 0q O = 0.65, and (b) 0^ o with y w = 0.475; for k = 0.02, 
L = 100 and two different values of y. 

of the CO surface, (ra), i.e. the relaxation time from the high CO coverage phase to the 
low-coverage phase, behaves in a similar way to (r p ). Figure ITUT a) clearly indicates that, to 
an even higher degree than in the poisoning process, (r^) is independent of y w . We choose 
y w = 0.57 for our calculations. Figure ITUTb) indicates that, as expected, (r d ) increases as 
9q decreases, but there is a range of values of the cut-off where the dependence is relatively 
small. In the following we calculate (rd) with 8q Q = 0.45, which lies in this region. 
We define the quantity 

A = \y-y 2 \, (9) 

which measures how far the system is from the coexistence curve and depends on k and 
L through y 2 . In Figs. ^2 and El we present snapshot configurations obtained during the 
relaxation from the low to the high coverage phase and from the high to the low coverage 
phase, respectively, when A -1 is small, i.e., far from the transition. In Figs. 1131 and 1141 we 
show snapshot configurations for a much larger value of A -1 , i.e., close to the transition. 
The difference in the decay mechanisms is evident from the figures. Figures ITT1 and [T2l clearly 
suggest that when A -1 is small, the system decays by nucleation and growth of multiple 
droplets of the stable phase. In contrast, Figs. and El show that when A -1 is large, 
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FIG. 10: (rd) as a function of (a) y w with 0q O 
L = 100 and two different values of y. 



0.45, and (b) 6* co with y w = 0.57; for k = 0.02, 



the system decays by nucleation and growth of a single droplet of the stable phase, which 
eventually takes over the entire system. The probability distributions of r p and rj, shown 
in Figs. HS1 and respectively, also indicate clear differences between the statistics of the 
decay times near and far from the coexistence line. Far from coexistence (A -1 small), the 
decay times follow an approximately Gaussian distribution. In contrast, near the transition 
(A -1 large), the distribution is approximately exponential. 

The statistics of the metastable lifetimes in the model studied here are strikingly similar 
to those found in Hamiltonian systems that decay toward thermodynamic equilibrium from 
a metastable phase associated with an egmlibrium phase transition. Well-studied examples 



are metastable decay in kinetic Ising (24 , 



25 



26 



27| and lattice-gas models [28, 



29( with such 



applications as magnetism switching and submonolayer adsorption. In the present paper we 
will for simplicity refer to this latter case as the "Hamiltonian" case, thus emphasizing the 
lack of a Hamiltonian and the consequent lack of a concept of thermodynamic equilibrium 
for the system studied in the present paper. The metastable decay in such a Hamiltonian 
system occurs via nucleation and subsequent growth of "droplets" inside which the order 
parameter is close to its equilibrium value, and it is well described by the classic KJMA 
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FIG. 11: Snapshot configurations obtained at different (unevenly spaced) times, after a sudden 
change of y from y w = 0.45 to y = 0.538 > y 2 (k), i.e., A^ 1 « 200. For k = 0.02 and L = 100. 
Here and in Figs. ^JE3 and^l the black dots represent lattice sites occupied by CO, while both 
adsorbed O atoms and empty lattice sites are white. 




FIG. 12: Snapshot configurations obtained at different (unevenly spaced) times, after a sudden 
change of y from y = 0.65 = y w to y = 0.5232 < y 2 (k), i.e., A -1 » 100. For k = 0.02 and L = 100. 

theory of phase transformation Q, The basic assumption of this theory is that 

droplets of the stable phase nucleate in a Poisson process £it £l r£ite / per unit volume. After 
nucleation, the droplet radius is assumed to grow with a constant speed, v. If the first 
droplet to nucleate grows fast enough to fill the system before another is likely to nucleate, 
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FIG. 13: Snapshot configurations obtained at different (unevenly spaced) times, after a sudden 
change of y from y = 0.45 = y w to y = 0.5338 > y 2 {k), i.e., A" 1 « 2000. For k = 0.02 and 
L = 100. 




FIG. 14: Snapshot configurations obtained at different (unevenly spaced) times, after a sudden 
change of y from y = 0.65 = y w to y = 0.5315 < y2(k), i.e., A -1 rs 500. For = 0.02 and L = 100. 

then it completes the phase transformation by itself - a process known as single-droplet 
(SD) decay. However, if the growth is slow, so that many droplets can nucleate within 
the time it would take a single droplet to fill the system, the phase transformation will 
proceed via a large number of droplets that nucleate and grow in parallel - a process known 
as multidroplet (MD) decay. This simple observation can be turned into a formal scaling 
argument by constructing the characteristic length R = (v/I) 1 ^ 3 , which is a measure of the 
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FIG. 15: Plot of distributions of r p for k = 0.02 and L = 100, (a) A" 1 « 200 and (b) A" 1 « 2000. 
Note the very different time scales in (a) and (b). 
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FIG. 16: Plot of distributions of r d for k = 0.02 and L = 100, (a) A^ 1 w 100 and (b) A" 1 w 500. 
Note the very different time scales in (a) and (b). 
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average size to which a droplet will grow before it touches another droplet. (Results are 



here reviewed only for d = 2. Results for general d can be found in, e.g., Refs. |2j, [26] and 
references therein.) If R S> L, the system is in the SD regime, and the metastable lifetime 
is simply the average time between nucleation events, rgo « 1/(IL 2 ). Since the nucleation 
events constitute a Poisson process, tsd follows an exponential distribution, similar to the 
ones shown in Figs. IToT b) and ITBT b) for the system discussed here. If, on the other hand, 
Ro <C L, then the system is in the MD regime, and the metastable lifetime is obtained as 
t M d = Ro/v = l/(f 2 /) 1//3 . Since the metastable decay in this case consists of a large number 
of droplets that nucleate and grow independently, tmd follows a Gaussian distribution with a 
standard deviation proportional to Rq/L, similar to the ones shown in Figs. IToT a) and lToT a) 
for the system discussed here. 

The SD regime with its exponential lifetime distribution is a subregion of a broader 
stochastic regime, while the MD regime with its narrow Gaussian distribution is part of a 
broader deterministic regime. The relative standard deviation of the lifetimes is defined by 

vV> - (r) 2 ( Ro/L < 1 in MD regime 
r = j- « < . (10) 

v / I 1 in SD regime 

The limit between the SD and MD regimes is called the Dynamic Spinodal (DSP) and 
corresponds to Ro ~ L. It is, however, easier to estimate it as given by the values of I and 
v that yield r = 1/2 Q. 

So far, the KJMA results discussed do not require a specific dependence of the nucleation 
rate I and growth velocity v on the macroscopic control parameters, which for metastable 
decay in Hamiltonian systems are the applied magnetic field H (or chemical potential or 
supersaturation for lattice-gas models) and the temperature T. In the present model the 
analogous quantities should be the distance from coexistence A and the desorption rate 
k. Hamiltonian systems are described by a free energy, and standard arguments of droplet 
theory show that for not too strong fields, / ~ exp [— c(T)/(T\H\)], where c(T) is well ap- 
proximated as proportional to the equilibrium interface tension between the metastable and 
equilibrium phases. For weak fields, v oc \H\ - an effect that to a reasonable approximation 
can be ignored compared to the exponential dependence on 1/\H\ in /. 

The present nonequilibrium system has no Hamiltonian and so no free-energy function. 
However, let us for the moment postulate that /(A, k) ~ exp(— c(k)/A) for reasonably small 
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FIG. 17: Log-linear plot of (a) (r p ) vs A -1 and (b) (rd) vs A" 1 , shown for fc = 0.01 and different 
values of L. The solid inverted triangles indicate AT,^ for each value of L. Approximate data 
collapse is seen for A" 1 < A-^. 



A, and that v(A, k) depends comparatively weakly on its parameters so that it can be taken 
as approximately constant. Following the method of data analysis introduced in Ref. 
we then expect logarithmic plots of (r p ) and (r^) (and of r in the MD regime) vs 1/A to be 
approximately linear. Furthermore, the general KJMA arguments given above indicate that 
the lifetimes should be independent of L in the MD regime and oc L~ 2 in the SD regime. 

Figure IT7I shows log-linear plots of (r p ) and (r^) vs A -1 for k = 0.01 and different sizes. 
Similar plots were also obtained for k = 0.02 (not shown). The plots clearly indicate that 
there is a regime, corresponding to small A~ x , where (r p ) is independent of L. Figure fTHl 
strongly indicates that when A -1 is large, the decay times are inversely proportional to 
1/L 2 . Only the data for L = 40 do not seem to follow this dependence. We believe it is 
possible that L = 40 becomes smaller than the critical droplet size for large A -1 (incipient 
"coexistence regime," see Ref. Q). 

To further explore the applicability of the KJMA theory and our postulate to the model, 
we next determine the dynamic spinodal, Adsp, that separates the stochastic and the de- 
terministic regimes. We calculate the relative standard deviation r of Eq. (jlOj) . which is 
shown on a logarithmic scale vs A -1 in Fig. EH where error bars are estimated by standard 
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FIG. 18: Log- linear plot of (a) L 2 (t p ) vs A -1 and (b) L 2 {ta) vs A -1 , shown for k = 0.01 and 
different values of L. The solid inverted triangles indicate A"^ for each value of L. Approximate 
data collapse is seen for A -1 > AT^. 



error-propagation methods as 



Or 



1 + 



n — 1 



n 



(11) 



As can be seen, r crosses over from the approximately linear behavior expected from our 
postulate for small A -1 to r ~ 1 for larger A -1 . We take as our estimate for Adsp the value 
Ax/2 for which r = 0.5. This crossover is determined for each value of L from the crossing 
of a weighted least-squares fit to In r in the linear region of Fig. with the horizontal line 
r = 0.5. The resulting estimates are shown in Fig. Ellas l/A^ vs L on a logarithmic scale, 
with error bars estimated from those in Fig.EUby standard error-propagation methods. The 
numerical results are consistent with the analytical prediction based on our postulate, 



^DSP 



1/hiL. 



(12) 



The asymptotic L dependence of the lifetime at the DSP, analogously given by 



(r) oc L/Adsp ~ L (InL) 



(13) 



is illustrated in Fig. For each L this lifetime was obtained by interpolation between the 
two closest field values bracketing Aj/2, for which simulations had been performed. The 
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FIG. 19: The relative standard deviation r of (a) (r p ) and (b) (rj), shown on a logarithmic scale 
vs A -1 for k = 0.01. The behavior of r crosses over from the approximate straight line expected 
from our postulate in the deterministic regime to r ~ 1 in the stochastic regime. The solid lines 
are weighted least-squares fits to the data in the linear region. 




40 60 100 150 200 




40 60 100 150 200 

L 



FIG. 20: The estimates l/A^ for 1/Adsp for (a) poisoning and (b) decontamination as obtained 
from Fig. EJand analogous data, shown vs L on a logarithmic scale. The solid lines are weighted 
least-squares fits excluding the points corresponding to L = 40. 
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FIG. 21: (a) (r p )/L at l/A^/g and (b) {t^)/L at I/A1/2, shown vs L on a logarithmic scale. The 
solid straight lines are weighted least-squares fits. Excluding the points corresponding to L = 40, 
the figures supports the asymptotic behavior (r p ) oc LlnL at the DSP. 



uncertainty in the resulting estimate was obtained by standard error propagation, taking 
the standard deviation in (r) from Eq. (jl(J|) with r = 1/2, and (r) from Fig. [*P71 

While a direct analogue of the surface tension does not exist in the present system, the 
results described above strongly suggest that it obeys a decay mechanism very similar to 
the one described by the standard KJMA theory of phase transformation by nucleation and 
growth, which predicts well-defined single-droplet and multidroplet regimes. A significant 
difference between our Fig. El and analogous figures showing the metastable lifetime for 
a Hamiltonian system vs inverse field or supersaturation (see, e.g, Fig. 2 of Ref. 24] and 
Fig. 2 of Ref. |27[), is that we here see no marked change in the slope of the curves at the 
DPT. One possible explanation is that the "effective surface tension" in the present case 
may decrease substantially with increasing A, in contrast to the situation in Hamiltonian 
systems. 

The decay times increase as y approaches the transition line 2/2 (k), however their behaviors 
depend on the direction of approach to the transition value, as can be seen in Fig. E21 I n a 



previous work we have shown how this asymmetry can be exploited to enhance the catalytic 



23 




0.5 0.51 0.52 0.53 0.54 0.55 0.56 

y 



FIG. 22: Decay times as functions of y when the system evolves toward the low CO coverage region 
(0q O = 0.45, y w = 0.55, left-pointing triangles) and when it evolves toward the high CO coverage 
region (Oqq = 0.65, y w = 0.475, right-pointing triangles) for k = 0.01 and k = 0.04, and L = 100. 
As discussed in the text, the divergence is exponential in l/\y — y2(k,L)\ = 1/A. 

activity by subjecting the system to periodic variation of the external pressure with periods 
related to the decay times in each direction [lo| . 

The asymmetry between the decay times when the system evolves to the low CO coverage 
phase, (r d ), and the decay time toward the high CO coverage phase, (r p ), becomes more 
evident as the desorption parameter k decreases, as can also be seen in Fig. E21 The extreme 
case occurs when k = 0, where (ra) = oo, independently of the value of y, since in the high 
CO coverage phase the surface then becomes irreversibly poisoned with CO, whereas (r p ) 
remains finite and dependent on y. 

In Fig. E2]it can be seen that (ra) and (r p ) appear to diverge at the same value of y. At 
this point the system spends, on average, the same amount of time in the low and high CO 
coverage phases, and it agrees with the value for the coexistence point y2, calculated from 
the order-parameter distribution in Sec. IIII Al Since at this point the average decay time 
is of the order of 10 4 MCSS, during the observation time of 10 6 MCSS it is very probable 
that several transitions occur between the two phases. These transitions are observed in 
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Fig. where it can be seen that the probability distribution, P(8co) is bimodal and quite 
symmetric, indicating that (ra) ~ (r p ) along the coexistence curve for finite L. 

IV. CONCLUSIONS 

We have investigated by kinetic Monte Carlo simulation the dynamical behavior of a 
ZGB model with desorption near the coexistence curve between the active and the CO 
poisoned nonequilibrium phases. We perform an extensive finite-size scaling analysis of the 
fluctuations and of the fourth-order reduced cumulant of the CO coverage, which plays the 
role of an order parameter. Our results strongly indicate, as also previously suggested by 
others, that the system undergoes a first-order nonequilibrium phase transition between the 
active and the CO poisoned phases. The coexistence curve terminates at a critical value of 
the desorption rate. We also calculated several points on the coexistence curve. 

Next we calculated the system-size dependence of the decay times of the metastable 
phases when the system is driven into the CO poisoned phase from the active phase, and 
vice versa. We found that near the coexistence curve the decay times are inversely propor- 
tional to 1/L 2 , and the decay mechanism consists of the nucleation and growth of a single 
supercritical droplet of the stable phase. In contrast, far from the coexistence curve, the 
decay times are independent of the system size, and the decay proceeds by random nucle- 
ation of many droplets of the stable phase that grow independently and coalesce. These 
regimes are separated by a dynamic spinodal that vanishes logarithmically with system 
size. These results strongly suggest that our nonequilibrium, non-Hamiltonian system fol- 
lows a decay mechanism very similar to the one described by the classic KJMA theory of 
phase transformation by nucleation and growth near a first-order equilibrium phase tran- 
sition, which predicts well-defined single-droplet and multidroplet regimes. In the present 
far-from-equilibrium system, the desorption parameter and the distance to the coexistence 
point play the roles of the temperature and the external field or supersaturation, respec- 
tively. Very recently, indications of KJMA behavior have also been observed in another 
non-Hamiltonian, nonequilibrium system: an ecological model of invasion by exotic species 
jsflj ] . We find quite exciting the strong similarity between the dynamics of metastable decay 
in far-from-equilibrium, non-Hamiltonian systems of applied importance, and the well-known 
behavior in systems that can be described by a Hamiltonian. 
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